library(tidyverse) # tidy universe
library(kableExtra) # html-table
library(patchwork) # combine plots
library(lmerTest) # mixed model
library(car) # ANOVA
library(performance) # model performanceConfounder age (Igf & animal welfare)
R Libraries
set.seed(1989)Data
Read data
Data processing done in file “igf-biomarker-summary-stats.qmd”.
load("./data/data-processed.RData")Statistical model: general design
With a simple model structure we solely test for a confounder effect of insemination age. We include a random animal effect and also the nested structure. Each sow gave birth 1, 2 or 3 times, so these events are nested within sows. We fit a linear mixed model with the lmerTest package in R.
contr = lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 10000000),
calc.derivs = FALSE)Cortisol serum
Model
mod.cort1 <- lmerTest::lmer(log(value) ~
insem.age +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "cort.ser") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)summary(mod.cort1)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "cort.ser")),
value)
Control: contr
REML criterion at convergence: 25.8
Scaled residuals:
Min 1Q Median 3Q Max
-1.88991 -0.48347 -0.01562 0.55338 1.54576
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 0.01358 0.1166
sow (Intercept) 0.02233 0.1494
Residual 0.05136 0.2266
Number of obs: 39, groups: litter.no:sow, 20; sow, 14
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 3.7595871 0.2312843 8.5866965 16.255 9.53e-08 ***
insem.age -0.0011780 0.0006513 8.0024430 -1.809 0.108
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
insem.age -0.965round(drop1(mod.cort1, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ insem.age + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age 0.168 0.168 1 8.002 3.271 0.108car::Anova(mod.cort1,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
insem.age 3.2711 1 0.07051 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.cort1,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
insem.age 2.7669 1 9.0492 0.1304Model diagnostics
performance::check_model(mod.cort1)
Plot
plot.cort1 <- dat.l |>
dplyr::filter(parameter == "cort.ser") |>
droplevels() |>
drop_na(value) |>
# make plot
mutate(jit = jitter(as.numeric(insem.age), 10)) |>
ggplot(aes(y = value)) +
geom_point(aes(x = jit,
col = husbandry,
shape = time),
size = 3) +
scale_colour_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 60),
breaks = seq(0, 60, 20),
minor_breaks = seq(0, 60, by = 2) ) +
labs(x = "Insemination age of the sows [days]",
y = "Cortisol (serum) [ng/ml]",
col = "Husbandry",
shape = "Time") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
col = guide_legend(override.aes = list(size = 5)))plot.cort1
Cortisol saliva
Model
mod.cort2 <- lmerTest::lmer(log(value) ~
insem.age +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "cort.sal") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.cort2)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "cort.sal")),
value)
Control: contr
REML criterion at convergence: 115.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.18019 -0.60262 0.02342 0.75178 1.69857
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 5.079e-17 7.127e-09
sow (Intercept) 0.000e+00 0.000e+00
Residual 8.584e-01 9.265e-01
Number of obs: 39, groups: litter.no:sow, 20; sow, 14
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.8128656 0.6986073 37.0000000 2.595 0.0135 *
insem.age -0.0004269 0.0019688 37.0000000 -0.217 0.8295
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
insem.age -0.977
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.cort2, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ insem.age + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age 0.04 0.04 1 37 0.047 0.83car::Anova(mod.cort2,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
insem.age 0.047 1 0.8284car::Anova(mod.cort2,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
insem.age 0.0402 1 11.608 0.8446Model diagnostics
performance::check_model(mod.cort2)
Plot
plot.cort2 <- dat.l |>
dplyr::filter(parameter == "cort.sal") |>
droplevels() |>
drop_na(value) |>
# make plot
mutate(jit = jitter(as.numeric(insem.age), 10)) |>
ggplot(aes(y = value)) +
geom_point(aes(x = jit,
col = husbandry,
shape = time),
size = 3) +
scale_colour_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 30),
breaks = seq(0, 30, 10),
minor_breaks = seq(0, 30, by = 2) ) +
labs(x = "Insemination age of the sows [days]",
y = "Cortisol (saliva) [ng/ml]",
col = "Husbandry",
shape = "Time") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
col = guide_legend(override.aes = list(size = 5)))plot.cort2
IGF bioactivity
Model
mod.bioact <- lmerTest::lmer(log(value) ~
insem.age +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "bioact.ser") |>
drop_na(value) |>
dplyr::filter(time != "30dpc") |>
droplevels(),
REML = TRUE,
control = contr)summary(mod.bioact)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
Data: droplevels(dplyr::filter(drop_na(dplyr::filter(dat.l, parameter ==
"bioact.ser"), value), time != "30dpc"))
Control: contr
REML criterion at convergence: 157.1
Scaled residuals:
Min 1Q Median 3Q Max
-2.41286 -0.33591 0.05923 0.42814 2.23561
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 0.2888 0.5374
sow (Intercept) 0.1375 0.3708
Residual 0.1364 0.3694
Number of obs: 80, groups: litter.no:sow, 40; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.391861 0.337537 33.159135 15.974 <2e-16 ***
insem.age -0.001524 0.000716 26.031980 -2.128 0.0429 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
insem.age -0.915round(drop1(mod.bioact, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ insem.age + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age 0.618 0.618 1 26.032 4.529 0.043 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.bioact,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
insem.age 4.5294 1 0.03332 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.bioact,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
insem.age 4.4425 1 27.39 0.04434 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics
performance::check_model(mod.bioact)
Plot
plot.bioact <- dat.l |>
dplyr::filter(parameter == "bioact.ser") |>
drop_na(value) |>
dplyr::filter(time != "30dpc") |>
droplevels() |>
# make plot
mutate(jit = jitter(as.numeric(insem.age), 10)) |>
ggplot(aes(y = value)) +
geom_point(aes(x = jit,
col = husbandry,
shape = time),
size = 3) +
scale_colour_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 600),
breaks = seq(0, 600, 200),
minor_breaks = seq(0, 600, by = 50) ) +
labs(x = "Insemination age of the sows [days]",
y = "IGF bioactivity [ng/ml]",
col = "Husbandry",
shape = "Time") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
col = guide_legend(override.aes = list(size = 5)))plot.bioact
IGF1 (serum)
Model
mod.igf1 <- lmerTest::lmer(log(value) ~
insem.age +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "igf1.ser") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.igf1)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igf1.ser")),
value)
Control: contr
REML criterion at convergence: 102.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.6324 -0.6741 0.1403 0.6757 1.5000
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 0.00441 0.0664
sow (Intercept) 0.00000 0.0000
Residual 0.45051 0.6712
Number of obs: 44, groups: litter.no:sow, 22; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 3.639462 0.471424 20.000005 7.720 2.01e-07 ***
insem.age 0.003985 0.001361 20.000005 2.928 0.00831 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
insem.age -0.976
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.igf1, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ insem.age + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age 3.863 3.863 1 20 8.574 0.008 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.igf1,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
insem.age 8.5744 1 0.003409 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.igf1,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
insem.age 7.4073 1 13.461 0.01702 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics
performance::check_model(mod.igf1)
Plot
plot.igf1 <- dat.l |>
dplyr::filter(parameter == "igf1.ser") |>
droplevels() |>
drop_na(value) |>
# make plot
mutate(jit = jitter(as.numeric(insem.age), 10)) |>
ggplot(aes(y = value)) +
geom_point(aes(x = jit,
col = husbandry,
shape = time),
size = 3) +
scale_colour_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 500),
breaks = seq(0, 500, 100),
minor_breaks = seq(0, 500, by = 20) ) +
labs(x = "Insemination age of the sows [days]",
y = "IGF1 (serum) [ng/ml]",
col = "Husbandry",
shape = "Time") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
col = guide_legend(override.aes = list(size = 5)))plot.igf1
IGF2 (serum)
Model
mod.igf2 <- lmerTest::lmer(log(value) ~
insem.age +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "igf2.ser") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.igf2)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igf2.ser")),
value)
Control: contr
REML criterion at convergence: 61.5
Scaled residuals:
Min 1Q Median 3Q Max
-3.2131 -0.8285 0.2732 0.6983 1.3886
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 0.0000 0.0000
sow (Intercept) 0.0000 0.0000
Residual 0.1937 0.4402
Number of obs: 40, groups: litter.no:sow, 20; sow, 14
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 4.451e+00 3.230e-01 3.800e+01 13.778 2.33e-16 ***
insem.age -7.232e-05 9.342e-04 3.800e+01 -0.077 0.939
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
insem.age -0.977
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.igf2, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ insem.age + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age 0.001 0.001 1 38 0.006 0.939car::Anova(mod.igf2,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
insem.age 0.006 1 0.9383car::Anova(mod.igf2,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
insem.age 0.005 1 12.7 0.9448Model diagnostics
performance::check_model(mod.igf2)
Plot
plot.igf2 <- dat.l |>
dplyr::filter(parameter == "igf2.ser") |>
droplevels() |>
drop_na(value) |>
# make plot
mutate(jit = jitter(as.numeric(insem.age), 10)) |>
ggplot(aes(y = value)) +
geom_point(aes(x = jit,
col = husbandry,
shape = time),
size = 3) +
scale_colour_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 200),
breaks = seq(0, 200, 50),
minor_breaks = seq(0, 200, by = 10) ) +
labs(x = "Insemination age of the sows [days]",
y = "IGF2 (serum) [ng/ml]",
col = "Husbandry",
shape = "Time") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
col = guide_legend(override.aes = list(size = 5)))plot.igf2
IGFBP2 (serum)
Model
mod.igfbp2 <- lmerTest::lmer(log(value) ~
insem.age +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "igfbp2.ser") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)summary(mod.igfbp2)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igfbp2.ser")),
value)
Control: contr
REML criterion at convergence: 72.2
Scaled residuals:
Min 1Q Median 3Q Max
-1.89723 -0.49844 0.00183 0.48081 2.04124
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 0.06597 0.2568
sow (Intercept) 0.07302 0.2702
Residual 0.05823 0.2413
Number of obs: 72, groups: litter.no:sow, 36; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.946e+00 1.871e-01 3.105e+01 31.773 <2e-16 ***
insem.age 8.381e-04 4.006e-04 2.365e+01 2.092 0.0473 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
insem.age -0.881round(drop1(mod.igfbp2, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ insem.age + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age 0.255 0.255 1 23.649 4.378 0.047 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.igfbp2,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
insem.age 4.3778 1 0.03641 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.igfbp2,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
insem.age 4.2439 1 24.209 0.0503 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics
performance::check_model(mod.igfbp2)
Plot
plot.igfbp2 <- dat.l |>
dplyr::filter(parameter == "igfbp2.ser") |>
droplevels() |>
drop_na(value) |>
# make plot
mutate(jit = jitter(as.numeric(insem.age), 10)) |>
ggplot(aes(y = value)) +
geom_point(aes(x = jit,
col = husbandry,
shape = time),
size = 3) +
scale_colour_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 1250),
breaks = seq(0, 1250, 250),
minor_breaks = seq(0, 1250, by = 100) ) +
labs(x = "Insemination age of the sows [days]",
y = "IGFBP2 (serum) [ng/ml]",
col = "Husbandry",
shape = "Time") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
col = guide_legend(override.aes = list(size = 5)))plot.igfbp2
IGFBP3 (serum)
Model
mod.igfbp3 <- lmerTest::lmer(log(value) ~
insem.age +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "igfbp3.ser") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.igfbp3)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igfbp3.ser")),
value)
Control: contr
REML criterion at convergence: 186.2
Scaled residuals:
Min 1Q Median 3Q Max
-1.4805 -0.6452 -0.1381 0.6411 1.7739
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 4.980e-01 7.057e-01
sow (Intercept) 2.263e-18 1.504e-09
Residual 3.266e-01 5.715e-01
Number of obs: 72, groups: litter.no:sow, 36; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 7.1029694 0.4448591 33.9999998 15.967 <2e-16 ***
insem.age -0.0001703 0.0009970 33.9999998 -0.171 0.865
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
insem.age -0.952
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.igfbp3, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ insem.age + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age 0.01 0.01 1 34 0.029 0.865car::Anova(mod.igfbp3,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
insem.age 0.0292 1 0.8644car::Anova(mod.igfbp3,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
insem.age 0.0281 1 27.386 0.8681Model diagnostics
performance::check_model(mod.igfbp3)
Plot
plot.igfbp3 <- dat.l |>
dplyr::filter(parameter == "igfbp3.ser") |>
droplevels() |>
drop_na(value) |>
# make plot
mutate(jit = jitter(as.numeric(insem.age), 10)) |>
ggplot(aes(y = value)) +
geom_point(aes(x = jit,
col = husbandry,
shape = time),
size = 3) +
scale_colour_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 6500),
breaks = seq(0, 6500, 2000),
minor_breaks = seq(0, 6500, by = 200) ) +
labs(x = "Insemination age of the sows [days]",
y = "IGFBP3 (serum) [ng/ml]",
col = "Husbandry",
shape = "Time") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
col = guide_legend(override.aes = list(size = 5)))plot.igfbp3
IGFBP2/IGFBP3 (serum)
Model
mod.igfbp23 <- lmerTest::lmer(log(value) ~
insem.age +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "igfbp23.ser") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.igfbp23)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igfbp23.ser")),
value)
Control: contr
REML criterion at convergence: 174.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.0444 -0.5111 0.1475 0.5307 1.6137
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 4.240e-01 6.512e-01
sow (Intercept) 7.539e-17 8.682e-09
Residual 2.736e-01 5.230e-01
Number of obs: 72, groups: litter.no:sow, 36; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -1.2353415 0.4096735 34.0000004 -3.015 0.00483 **
insem.age 0.0010469 0.0009182 34.0000004 1.140 0.26215
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
insem.age -0.952
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.igfbp23, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ insem.age + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age 0.356 0.356 1 34 1.3 0.262car::Anova(mod.igfbp23,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
insem.age 1.3001 1 0.2542car::Anova(mod.igfbp2_3,
test.statistic = "F",
type = 2)Error: Objekt 'mod.igfbp2_3' nicht gefundenModel diagnostics
performance::check_model(mod.igfbp23)
Plot
plot.igfbp23 <- dat.l |>
dplyr::filter(parameter == "igfbp23.ser") |>
droplevels() |>
drop_na(value) |>
# make plot
mutate(jit = jitter(as.numeric(insem.age), 10)) |>
ggplot(aes(y = value)) +
geom_point(aes(x = jit,
col = husbandry,
shape = time),
size = 3) +
scale_colour_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 2.5),
breaks = seq(0, 2.5, 0.5),
minor_breaks = seq(0, 2.5, by = 0.1) ) +
labs(x = "Insemination age of the sows [days]",
y = "Molar ratio IGFBP2/IGFBP3 (serum)",
col = "Husbandry",
shape = "Time") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
col = guide_legend(override.aes = list(size = 5)))plot.igfbp23
Proteolytic activity (serum)
Model
mod.prot <- lmerTest::lmer(log(value+1) ~
insem.age +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "proteolysis") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.prot)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value + 1) ~ insem.age + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "proteolysis")),
value)
Control: contr
REML criterion at convergence: 132.9
Scaled residuals:
Min 1Q Median 3Q Max
-1.5185 -1.2360 0.2984 0.8499 1.3739
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 0.000 0.000
sow (Intercept) 0.000 0.000
Residual 1.339 1.157
Number of obs: 39, groups: litter.no:sow, 20; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.205e+00 5.616e-01 3.700e+01 2.145 0.0386 *
insem.age 8.813e-04 1.406e-03 3.700e+01 0.627 0.5345
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
insem.age -0.944
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.prot, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value + 1) ~ insem.age + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age 0.526 0.526 1 37 0.393 0.535car::Anova(mod.prot,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value + 1)
Chisq Df Pr(>Chisq)
insem.age 0.3931 1 0.5307car::Anova(mod.prot,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value + 1)
F Df Df.res Pr(>F)
insem.age 0.3027 1 15.645 0.59Model diagnostics
performance::check_model(mod.prot)
Plot
plot.prot <- dat.l |>
dplyr::filter(parameter == "proteolysis") |>
droplevels() |>
drop_na(value) |>
# make plot
mutate(jit = jitter(as.numeric(insem.age), 10)) |>
ggplot(aes(y = value)) +
geom_point(aes(x = jit,
col = husbandry,
shape = time),
size = 3) +
scale_colour_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 21),
breaks = seq(0, 21, 5),
minor_breaks = seq(0, 21, by = 1) ) +
labs(x = "Insemination age of the sows [days]",
y = "Proteolytic activity (serum) [%]",
col = "Husbandry",
shape = "Time") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
col = guide_legend(override.aes = list(size = 5)))plot.prot
Stanniocalcin, STC1 (serum)
Model
mod.stc1 <- lmerTest::lmer(log(value) ~
insem.age +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "stc1") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)summary(mod.stc1)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "stc1")),
value)
Control: contr
REML criterion at convergence: 105.4
Scaled residuals:
Min 1Q Median 3Q Max
-1.59654 -0.47741 -0.00781 0.45021 1.56136
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 2.0946 1.447
sow (Intercept) 1.6464 1.283
Residual 0.0882 0.297
Number of obs: 38, groups: litter.no:sow, 19; sow, 14
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.264645 1.796520 5.694156 1.261 0.2567
insem.age 0.012553 0.005163 5.230989 2.431 0.0571 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
insem.age -0.963round(drop1(mod.stc1, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ insem.age + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age 0.521 0.521 1 5.231 5.911 0.057 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.stc1,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
insem.age 5.9111 1 0.01505 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.stc1,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
insem.age 4.8382 1 7.3349 0.06203 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics
performance::check_model(mod.stc1)
Plot
plot.stc1 <- dat.l |>
dplyr::filter(parameter == "stc1") |>
droplevels() |>
drop_na(value) |>
# make plot
mutate(jit = jitter(as.numeric(insem.age), 10)) |>
ggplot(aes(y = value)) +
geom_point(aes(x = jit,
col = husbandry,
shape = time),
size = 3) +
scale_colour_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 10000),
breaks = seq(0, 10000, 2000),
minor_breaks = seq(0, 10000, by = 400) ) +
labs(x = "Insemination age of the sows [days]",
y = "STC1 (serum) [ng/ml]",
col = "Husbandry",
shape = "Time") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
col = guide_legend(override.aes = list(size = 5)))plot.stc1
Calcium (serum)
Model
mod.calc <- lmerTest::lmer(log(value) ~
insem.age +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "calc") |>
dplyr::filter(time != "30dpc") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.calc)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ insem.age + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dplyr::filter(dat.l, parameter ==
"calc"), time != "30dpc")), value)
Control: contr
REML criterion at convergence: 11.7
Scaled residuals:
Min 1Q Median 3Q Max
-5.9541 -0.2379 0.0437 0.3503 1.7735
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 4.477e-19 6.691e-10
sow (Intercept) 3.601e-02 1.898e-01
Residual 3.934e-02 1.983e-01
Number of obs: 80, groups: litter.no:sow, 40; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.388e+00 9.118e-02 6.227e+01 26.190 <2e-16 ***
insem.age 2.126e-05 1.720e-04 6.730e+01 0.124 0.902
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
insem.age -0.804
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.calc, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ insem.age + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
insem.age 0.001 0.001 1 67.303 0.015 0.902car::Anova(mod.calc,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
insem.age 0.0153 1 0.9016car::Anova(mod.calc,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
insem.age 0.0151 1 25.64 0.9032Model diagnostics
performance::check_model(mod.calc)
Plot
plot.calc <- dat.l |>
dplyr::filter(parameter == "calc") |>
dplyr::filter(time != "30dpc") |>
droplevels() |>
drop_na(value) |>
# make plot
mutate(jit = jitter(as.numeric(insem.age), 10)) |>
ggplot(aes(y = value)) +
geom_point(aes(x = jit,
col = husbandry,
shape = time),
size = 3) +
scale_colour_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 21),
breaks = seq(0, 20, 5),
minor_breaks = seq(0, 20, by = 2) ) +
labs(x = "Insemination age of the sows [days]",
y = "Calcium (serum) [mg/dl]",
col = "Husbandry",
shape = "Time") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
col = guide_legend(override.aes = list(size = 5)))plot.calc
Combined plots: Figure
# Combine plots with a designated area for the legend
combined <- (plot.cort1 +
plot.cort2 +
plot.bioact +
plot.igf1 +
plot.igf2 +
plot.igfbp2 +
plot.igfbp3 +
plot.igfbp23 +
plot.prot +
plot.stc1 +
plot.calc +
guide_area()) +
plot_layout(guides = "collect") +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(face = "bold", size = 20),
legend.position = "right",
legend.direction = "vertical")combined
png("./plots/figure-age.png",
width = 600, height = 500, units = "mm",
pointsize = 10, res = 600)
combined
dev.off()png
2 How to cite R
“All analyses were performed using R Statistical Software (version 4.4.2; R Core Team 2024)”.
Reference: R Core Team (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
citation()To cite R in publications use:
R Core Team (2024). _R: A Language and Environment for Statistical
Computing_. R Foundation for Statistical Computing, Vienna, Austria.
<https://www.R-project.org/>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Manual{,
title = {R: A Language and Environment for Statistical Computing},
author = {{R Core Team}},
organization = {R Foundation for Statistical Computing},
address = {Vienna, Austria},
year = {2024},
url = {https://www.R-project.org/},
}
We have invested a lot of time and effort in creating R, please cite it
when using it for data analysis. See also 'citation("pkgname")' for
citing R packages.version$version.string[1] "R version 4.4.2 (2024-10-31 ucrt)"citation("performance")Um Paket 'performance' in Publikationen zu zitieren, nutzen Sie bitte:
Lüdecke et al., (2021). performance: An R Package for Assessment,
Comparison and Testing of Statistical Models. Journal of Open Source
Software, 6(60), 3139. https://doi.org/10.21105/joss.03139
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Article{,
title = {{performance}: An {R} Package for Assessment, Comparison and Testing of Statistical Models},
author = {Daniel Lüdecke and Mattan S. Ben-Shachar and Indrajeet Patil and Philip Waggoner and Dominique Makowski},
year = {2021},
journal = {Journal of Open Source Software},
volume = {6},
number = {60},
pages = {3139},
doi = {10.21105/joss.03139},
}Session Info
sessionInfo()R version 4.4.2 (2024-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)
Matrix products: default
locale:
[1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8
[3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
[5] LC_TIME=German_Germany.utf8
time zone: Europe/Berlin
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] performance_0.13.0 car_3.1-3 carData_3.0-5 lmerTest_3.1-3
[5] lme4_1.1-36 Matrix_1.7-1 patchwork_1.3.0 kableExtra_1.4.0
[9] lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[13] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
[17] ggplot2_3.5.1 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] gtable_0.3.6 bayestestR_0.15.1 xfun_0.50
[4] htmlwidgets_1.6.4 ggrepel_0.9.6 insight_1.0.1
[7] lattice_0.22-6 tzdb_0.4.0 numDeriv_2016.8-1.1
[10] vctrs_0.6.5 tools_4.4.2 Rdpack_2.6.2
[13] generics_0.1.3 datawizard_1.0.0 parallel_4.4.2
[16] pbkrtest_0.5.3 pkgconfig_2.0.3 lifecycle_1.0.4
[19] compiler_4.4.2 farver_2.1.2 munsell_0.5.1
[22] htmltools_0.5.8.1 yaml_2.3.10 Formula_1.2-5
[25] pillar_1.10.1 nloptr_2.1.1 MASS_7.3-61
[28] reformulas_0.4.0 boot_1.3-31 abind_1.4-8
[31] nlme_3.1-166 tidyselect_1.2.1 digest_0.6.37
[34] stringi_1.8.4 labeling_0.4.3 splines_4.4.2
[37] fastmap_1.2.0 grid_4.4.2 colorspace_2.1-1
[40] cli_3.6.3 magrittr_2.0.3 broom_1.0.7
[43] withr_3.0.2 backports_1.5.0 scales_1.3.0
[46] timechange_0.3.0 rmarkdown_2.29 hms_1.1.3
[49] evaluate_1.0.3 knitr_1.49 rbibutils_2.3
[52] viridisLite_0.4.2 mgcv_1.9-1 rlang_1.1.4
[55] Rcpp_1.0.13-1 glue_1.8.0 xml2_1.3.6
[58] see_0.9.0 svglite_2.1.3 rstudioapi_0.17.1
[61] minqa_1.2.8 jsonlite_1.8.9 R6_2.6.1
[64] systemfonts_1.1.0